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Abstract 

The thermal conductivity of a crystal is sensitive to the presence of surfaces and nanoscale defects. 
While this opens tremendous opportunities to tailor thermal conductivity, a true "phonon engi- 
neering" of nanocrystals for a specific electronic or thermoelectric application can only be achieved 
when the dependence of thermal conductivity on the defect density, size, and spatial population 
is understood and quantified. Unfortunately, experimental studies of effects of nanoscale defects 
are quite challenging. While molecular dynamics simulations are effective in calculating thermal 
conductivity, the defect density range that can be explored with feasible computing resources is 
unrealistically high. As a result, previous work has not generated a fully detailed understanding 
of the dependence of thermal conductivity on nanoscale defects. Using GaN as an example, we 
have combined physically-motivated analytical model and highly-converged large scale molecular 
dynamics simulations to study effects of defects on thermal conductivity. An analytical expression 
for thermal conductivity as a function of void density, size, and population has been derived and 
corroborated with the model, simulations, and experiments. 
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I. INTRODUCTION 



Ability to tailor thermal transport properties of materials has become increasingly im- 
portant in modern society. In one example, reducing thermal conductivity of thermoelectric 
materials makes it practical to generate electricity from waste heat [HE]- In another ex- 
ample, increasing thermal conductivity of the semiconductor materials enables an increase 
in the device density by improving heat dissipation [3j H]. It is well known that thermal 
conductivities of materials are extremely sensitive to defect populations [5T4TT]. Control of 
defects therefore provides a powerful means to minimize thermal conductivities for thermo- 
electric applications. The presence of defects, however, is a concern to the electronic devices 
especially in the nano-scale. This is because defects not only affect electronic properties, 
but also cause phonon scattering that reduces thermal conductivities. This can lead to 
abnormal heating in local regions around defects. If the heated region is the bottle-neck 
of a nanostructure, for example, a cross-section of a nanowire, catastrophic failure may be 
triggered [12]. A detailed understanding of thermal conductivity as a function of defects is 
thus essential. 

Given experimental evidence of significant defect incorporation p~3HTT| . studies of defects 
in GaN nanostructures are of particular interest. Specifically, recent experiments [18] suggest 
that the spatial distribution of defects in GaN nanowires is not uniform, and vacancies can 
combine to form nano-scale voids near surfaces. Understanding thermal conductivity as a 
function of void population using GaN as an example is therefore technologically beneficial 
since GaN is an important semiconductor that can be readily integrated with Si and is 
ideally suited for optoelectronic applications [19J. 

Various experiments have been performed to quantify thermal conductivity as a function 
of porosity of porous materials [201421] . However, experimental studies of nanoscale voids 
are much more challenging and a detailed knowledge of void effects has not been developed 
in experiments. Extensive theoretical work has also been performed to develop relationships 
between thermal conductivity and point defects [7J El 1251 - 130] . Most of these studies [7J El I25T 
[29] used an analytical approach initially developed by Callaway and Vonbaeyer [25] , Klemens 
[26] . Abeles [27], and Slack [28] et al. In this approach, a point defect phonon scattering 
relaxation time is expressed as a function of the atomic masses, and the strain caused by the 
defects. This relaxation time is then used in Debye model to calculate thermal conductivity. 
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While this approach has been successfully applied to explore thermal conductivity as a 
function of defect density, the defect scattering parameter in the model is often treated as a 
fitting parameter requiring experimental data to determine. Also, this model has not been 
used to distinguish different types of defects such as vacancies and voids (with different 
sizes). In particular, we note that the model was mainly used for point defects but not for 
bigger defects such as voids. 

By using a discrete lattice directly in a computational system where defects of different 
sizes and shapes can be precisely replicated on an atomistic scale, molecular dynamics (MD) 
enables the effects of voids on thermal conductivity to be predicted from a known interatomic 
potential without additional assumptions. MD therefore provides a powerful means to study 
thermal transport [30HiO] . In particular, MD has been used to study the effects of vacancy, 
interstitial, and antisite defects on thermal conductivity [301 1381H0] , and to determine the 
scattering of phonon wave packets by point defects [JT]. Using relatively small systems 
imposed by the computational cost, these previous MD studies were limited to very high 
defect density (site fraction 0.5% [3D], > 0.0125% [38J, > 0.39% [39J, > 0.075% [ID], and > 
0.016% |41j) even when a single defect is simulated. In contrast, a very high experimental 
vacancy concentration of 10 15 cm -3 only corresponds to a site fraction of 1.15 x 10~ 6 %. It 
is unclear how the MD results scale with the experimental data obtained at defect densities 
orders of magnitude smaller. 

Recently, we have developed and verified a scaling law |42j that can extend the size- 
dependent MD conductivity estimates to the dimensions of realizable devices. In particular, 
analytical expressions of thermal conductivity as a function of film thickness or wire radius 
are derived [43J . Such a scaling law provides an ideal means to evaluate the scaling of defect 
density from MD to experimental conditions (a discussion of the validity of the scaling law 
is given in Appendix [A]) . Here the scaling law of defects can be explored using three defect 
distributions as shown in Fig. [I] In Fig. [TJa), the defects are uniformly distributed in all 
three directions with equal spacing 5. In Fig. [lib), defects are closely spaced in the x— 
and z— directions with spacings of i x and i z respectively but the x — z defect arrays have 
a large spacing 5 in the y— direction. In Fig. [IJc), defects are closely spaced in the y— 
and z— directions with spacings of i y and i z respectively but the y — z defect arrays have a 
large spacing 5 in the x— direction. Increasing 5 can effectively reduce the void site fraction. 
Despite different distributions, the mechanism and the functional dependence of thermal 
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conductivity on 5 is similar in all three cases if the spacings in each defect array (i.e., l x , t y , 
l z ) remain constant. Compared with Fig. [ija), smaller systems can be used in Figs. []Jb) 
and[T^c) to capture the same defect density. Hence, configurations shown in Figs. Sib) and 
[TJc) are more suitable for the computationally extremely expensive MD simulations. 

The present work addresses the scaling of defect effects on thermal conductivity in the 
realistic defect density range using three integrated approaches: (a) develop a physically- 
motivated model and use it to derive an analytical expression of thermal conductivity as a 
function of void density, size, and spatial population; (b) verify the analytical model using 
massively parallel MD simulations in a case study that studies the effects of voids on the 
thermal conductivity of GaN wurtzite crystal; and (c) corroborate the analytical model 
and MD simulations using the extensive experimental thermal conductivity data of porous 
materials. The significance of the recent experimental observation that the nano-voids in 
some GaN nanowires appear mostly near the surfaces rather than in the bulk [18] will also 
be evaluated. 

II. MOLECULAR DYNAMICS SIMULATION METHODOLOGY 

Molecular dynamics simulations are used to both develop and verify the analytical model 
of defects. In particular, the "direct method" [33J 1121 S3] was used to calculate the [0001] 
thermal conductivity / resistivity of a defective wurtzite GaN crystal at a temperature of 
T = 300 K. Here the temperature was defined as T = 2ek/{3k), where is average kinetic 
energy per atom and k is Boltzmann constant. We caution that the 300 K temperature is at 
the lower bound of the estimated Debye temperature range of 350 - 600 K [HH16], and we 
did not perform quantum correction on temperature as it has been found to be of dubious 
benefit [47J. However, the use of the low temperature is satisfactory in the present work 
because our purpose is to explore the scaling of thermal conductivity with defect density 
rather than evaluate the accuracy of the MD method on thermal conductivity calculations 
at low temperatures. One particular reason to choose the relatively low temperature is to 
reduce intrinsic noise of the data j33J S21 S3] so that the scaling relationship can be more 
reliably identified. Previous MD work [331 S21 S3] has employed the Stillinger- Weber (SW) 
potential parameterized by Bere and Serra [Ml EH] to calculate the thermal conductivity of 
GaN bulk and nanostructured crystals. To compare with the previous results, we use the 
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FIG. 1: Three void distributions: (a) a bulk configuration with a uniform density; (b) a configura- 
tion with arrays of voids on the x-z plane; and (c) a configuration with arrays of voids on the y-z 
plane. 

same potential in the present study. This potential gives reasonable prediction on dispersion 
relations, vibrational density of states (DOS), and heat capacity for the bulk GaN system 
[33] . To ensure accuracy of our MD thermal conductivity data, all of our calculations are 
based upon highly converged temperature profiles averaged over at least 50 ns, which is 
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significantly longer than in similar studies. 

The computational system used for the simulations is shown in Fig. [2] The hexagonal 
wurtzite GaN crystal is aligned so that the x—, y—, and z— coordinates correspond respec- 
tively to [0001], [1100], and [1120] directions. The computational system has a length 2L 
in the x— direction, a thickness t in the y— direction, and a width W in the z— direction. 
According to the lattice constants of GaN and the chosen geometry, the smallest orthogonal 
cell has a dimension of a\ — c — 5.2000 A, a 2 = 2 ■ a ■ cos (7r/6) = 5.5252 A, and = a = 
3.1900 A in the x—, y—, and z— directions [23] • For convenience, the system dimensions 
are represented by the number of cells riz, and n 3 in the x—, y—, and z— directions. A 
constant heat flux method [321 |33j [35], H2} 03j I50H52"] is used to create a temperature gradient 
in the x— direction as shown by the color scheme in Fig. [2] (red and blue mean respectively 
the highest and the lowest temperatures). Void defects are created uniformly between heat 
source and sink by removing clusters of atoms. In particular, the void is equivalent to a 
volume of 12 atoms and a shape of hexagonal column along the x— direction as shown in 
the blow-up of Fig. [2] Such hexagonal voids are terminated with low index surfaces and are 
stable during simulations. 

Thermal transport (along x— direction) in both bulk and film configurations is considered. 
As shown in Fig. |2j the bulk means that the system thickness in the y— direction is infinite, 
and this is achieved by using a periodic boundary condition in the y— direction. The film 
means that the system thickness in the y— direction is finite, and this is achieved by using a 
free boundary condition in the y— direction. The periodic boundary conditions are always 
used in the x— and z— directions for both bulk and film configurations. It should be noted, 
however, that the phonon mean free path in the x— direction is limited by the finite distance 
between the heat source and sink despite the periodic boundary condition. 

III. ANALYTICAL MODEL OF THE THERMAL CONDUCTIVITY OF MATE- 
RIAL WITH VOIDS 

To ensure that our analytical model of defects is physically well-motivated, MD simula- 
tions are first performed to examine the effect of defects on the heat flux distribution. A film 
configuration composed of around 590,000 atoms with dimensions of n\ = 136 (2L ~ 707A), 
n,2 = 90 (t m 500 A), and n$ = 6 (W ~ 19 A) is used. Following the same methodology 
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FIG. 2: Geometry of "direct method" MD simulations of thermal conductivity. 

described previously [13] , a local heat flux is calculated as a long time average of the atomic 
contributions to the overall heat flux [32] for atoms in the rectangular box illustrated in Fig. 
^a). The box has 150 A x— dimension and about 41 A y— dimension, and extends all 
the z— dimension of computational cell. It is placed at an x— coordinate of about 150 A. 
By moving the box along the y— direction, heat flux is calculated as a function of y. Three 
simulations, corresponding respectively to: no defects, four voids equally spaced between 
the heat source and sink along the center plane of the film, and the same four voids along a 
plane near the y— surface of the film (precisely 2 unit cells from the surface), are conducted 
as shown in Fig. ga). The heat fluxes obtained from the three simulations are shown in 
Fig. gb). 

Fig. gb) indicates that for the sample with no defects, the heat flux inside the film 
remains relatively constant and the heat flux decreases only near the surfaces due to the 
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(a) system geometry and location of flux measurement 
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FIG. 3: Heat flux distributions in the cross-section of a thin film: (a) system geometry illustrating 
three void configurations (1: no voids; 2: voids near center; and 3: voids near one of the surfaces), 
box for calculating the local heat flux, and the x— coordinate where the box is placed and moved 
in the y— direction to calculate heat flux as a function of y— coordinate; and (b) heat flux as a 
function of y— coordinate for the three void configurations. The heat flux is normalized by the 
average flux value to emphasize the relative differences in flux across the film. 
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surface scattering effect, as reported in previous work [13]. From the graph, the scattering 
distance of the surface-affected regions can be estimated to be about 100 A. This is in 
agreement with previous molecular dynamics calculations of thermal conductivities of wires 
and films [121 S3] which indicated that at least for the GaN SW potential [A8\ H9] , the surface 
scattering affected heat flux is confined to a near-surface layer around 100 A deep. Fig. [3] 
also indicates that voids in the center of the film causes a reduction of heat flux due to 
the defect scattering effect. Apparently, the dimension of the affected region is also about 
100 A (or less) for voids of this size. The heat flux in regions away from the voids and 
the surfaces is nearly constant and comparable to the values of the film without a defect. 
When the voids are moved to near the surface, the void and surface scattering zones are 
consolidated so that only a combined reduction of heat flux near the surface is observed. 
The difference is apparent in the reduction at the surface with the near-surface void which 
is slightly more significant than that at the other surface. These findings are the basis for 
the analytical model of the heat flow to be developed in the following. 

Phonon density of states is also explored. The same system dimension as that in Fig. [3] 
is used with four voids equally spaced between heat source and sink along the center plane. 
Periodic boundary condition is used in the y— direction to eliminate the surface effects so 
that the void effects can be more clearly revealed. Phonon density of states is calculated for 
five selected regions shown in Fig. 121(a) using the real scale (red boxes), where region 1 to 
region 5 has increasing distance from the voids. From the calculated density of states shown 
in Fig. |4^b), it is apparent that there is no significant change in the acoustic density of states 
across the system and there is a subtle increase in the optical density of states. The increase 
in optical states near the void is plausible. However, a full understanding of the influence of 
voids on individual phonons and the contributions of phonons to the observed flux behavior 
would require analysis of local mean free paths that is beyond the present scope [53] . 

A. Introduction of Model Concepts 

According to Fig. j3j we assume that the mean effects of phonon scattering at surfaces 
and voids is confined to the 100 A ranges, at least for voids of comparable sizes to those 
simulated. Considering that the average spacing between point defects such as vacancies 
exceeds 1000 A even when the material contains a high defect concentration of 10 15 cm~ 3 , 



9 
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FIG. 4: The phonon density of states for a periodic system: (a) locations of 5 regions (red boxes) 
for calculations of density of states (from region 1 closest to the voids to region 5 furthest); and 
(b) the density of states calculated for the 5 locations shown in (a). Here the values of density 
of states have error bars of near the width of the variation of the acoustic modes (0-45 meV) and 
are omitted for clarity. The density of acoustic modes appears to be equivalent but the density of 
optical modes increases slightly near the void. 
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we can assume that defects are dispersed (i.e., isolated and independent). For convenience, 
we further assume that defects are uniformly distributed, which is roughly true for the most 
homogeneous experimental material, but is exact for simulations with periodic boundary 
conditions where defects are replicated uniformly in space. On the other hand, when the 
defect scattering effect is fully confined to a relatively small region, the overall thermal 
conductivity becomes independent on the particular distribution of sparse defects. This is 
examined in Fig. [5] using heat transport through a prism as an example. For sparse defects, 
we can embed each defect (circle) in a section (dark box) with a sufficiently large dimension 
of d x . Because defect effects are confined to a relatively small region, we can assume that 
only the thermal conductivity inside the section changes to whereas the conductivity 
outside the section remains to be the defect-free matrix value K m . It is plausible that, for 
the two cases where the three defects are either uniformly distributed over the prism length 
L, Fig. |5](a), or concentrated on the left-hand side, Fig. ^b), the overall thermal resistivities 
(inverse of thermal conductivity) of the two prisms are the same, both equal to a length 
weighted average of resistivities of defective sectors and matrix [12J H3] as: 3d x /L ■ n^ 1 + 
(L — 3d x ) I L ■ k" 1 . Note that our model will be accurate when these assumptions are valid. 
The model may still capture well the scaling law even for materials with phonon mean free 
path longer than defect spacing. Model evaluation should be performed by comparing model 
predictions with either experiments or direct simulations (such as MD simulations as will 
be shown in the following for the GaN case). 

When the defect distribution is uniform, the material can be viewed as composed of 
equivalent small volumes 5V — 6 X ■ 8 y ■ S z each containing a defect. Because these volumes 
are assumed to be identical and independent of each other, the overall thermal conductivity 
of the material equals the thermal conductivity of each individual volume. As a result, 
only one representative volume needs to be addressed. Here, we consider heat conduction 
in the x dimension of the volume 5V = 5 X ■ 5 y ■ 5 Z containing a defect in its center as 
marked by a spherical ball in Fig. |6](a). The presence of this defect changes the thermal 
conductivity of its surroundings defined by a sub- volume dV = d x -d y -d z . Because the defect 
concentration is low, we can always choose large values of d x , d y , and d z (comparable to the 
defect scattering distance shown in Fig. [3]). Hence, the material outside the dV = d x ■ d y ■ d z 
volume can be viewed as far away from the defect and therefore its thermal conductivity 
remains to be the value of the (defect-free) matrix material, n m . Despite the non-uniform 
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(a) uniform defect distribution 

s dx ^ heat flux ► 



K m 








K m 


Q K d 


K m 


< 












> 



L 

(b) nonuniform defect distribution 



o 


o 


^9 





FIG. 5: Comparison of two defect distributions: (a) uniform and (b) nonuniform defect distribu- 
tions. 

thermal properties on a fine scale, the smaller scattering volume dV exhibits an apparent 
overall thermal conductivity k^. Since is essentially the coarse-grained, average thermal 
conductivity of the volume d x -d y -d z , it therefore depends strongly on d x , d y , and d z ; however, 
can be viewed as a constant for a given defect size and type once d x , d y , and d z are given. 

Using the model shown in Fig. [6^a), we define two geometric parameters a = 
(d y ■ d z ) f(8 y ■ S z ) and (3 = d x /5 x . Considering that S x , S y , and 5 Z are the overall mate- 
rial dimensions divided by numbers of defects in the three coordinate directions, and d x , d y , 
and d z are given constants physically representing the dimension around the defect where 
scattering is significant, parameters a and j3 prescribe the relative defect densities (areal and 
lineal fractions) in the cross-sectional and axial dimensions respectively. In particular, a and 
(3 can be termed defect "scattering" densities because a = 1 and /3 — 1 do not represent 
100% defect volume fraction (or site fraction), but rather indicate that the ratio between 
defect scattering volume and total volume equals one. We can also define a defect volume 
scattering density p = a - (3. Note that for a given defect size and type, the defect scattering 
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FIG. 6: Heat conduction through a volume S x ■ S y ■ S z of material containing one defect: (a) average 
total volume (Sx ■ 5y ■ Sz) containing one defect; (b) prism (Sx ■ dy ■ dz) containing one defect; and 
(c) slab (dx ■ 5y ■ Sz) containing one defect. 
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density p is proportional to the defect volume fraction £, £ = / • p, with the scaling factor / 
defined as 

Q P ■ P ■ P 

; = n~ 8 = (4 + Q ■ (£ y + 4) ■ (4 + is) 

where Q is physical volume of defect, fl s represents a scattering volume, P x , P y , and i z 
are physical dimensions of the defect in the x—, y—, and z— directions, and £ s is a mean 
scattering distance for all phonons contributing to the heat flux. Eq. indicates that / 
depends on defect size. As will be clear below, this is an important concept illustrating why 
thermal conductivity is sensitive to defect sizes. 

Our distinction of two defect densities a and f3 may not appear critical for uniform defect 
distributions. These two separate measures of defect density, however, are useful when the 
defect distribution is not uniform. For instance, consider the case of heat conduction through 
a long nanowire that is almost defect-free along its axial direction (/3 « 0) but is severely 
damaged or even completely fractured on one cross-section (a ~ 1). Such a nanowire is 
expected to have a low thermal conductivity. On the other hand, if the same amount of 
defects is uniformly distributed along the axis of a long nanowire, defect densities become 
low in both cross-sectional and axial dimensions (a ~ 0, /? ~ 0). Such a nanowire is expected 
to have a thermal conductivity close to that of the matrix, K m . Clearly, the two parameters 
a and (3 can be used to distinguish the two cases. 



B. Derivation of Thermal Conductivity as a Function of Defect Density, Size, and 
Distribution 

Analytical expression can be derived for thermal conductivity of a material as a function 
of defect density, size, and distribution using the coarse-grained conductivity concept just 
introduced. The free parameters of this model are the matrix and defect volume conductivi- 
ties, K m and Kd, and the chosen defect volume size d x , d y , and d z . First, the material volume 
8 x -8 y - 5 Z is divided along the x direction into a central rectangular prism with a cross-section 
area d y ■ d z that contains the defect and the surrounding material that is far away from the 
defect, as shown in the left frame of Fig. ^b). The apparent overall thermal conductivity of 
the central prism can be assumed to be Kd y d z , whereas that of the remaining material is K m . 
Because the heat conducts through the two parts of the material in parallel (on average), 
the thermal conductivity of the entire material (matrix + defect) can be calculated as an 
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area weighted average [321 133] : 



K m+d = T J- ' K dyd* H x~i~ K m = a * «A,d, + (1 " tt) • K m (2) 



In Eq. (|2j), w^d, can be expanded by observing that the central prism is composed of a 
central section with a length d x and a conductivity K d and two end sections with a total 
length 5 X — d x and a conductivity K m as shown in the right frame of Fig. (6^b) . Because heat 
conducts through the three sections in serial, the overall thermal resistivity is calculated as 
a length weighted average [321 S3] : 

1 d x 1 fi x - d x 1 1-/3 . s 



K dyd z fix K d fix K m ^d K 

Substituting Eq. ([3]) into Eq. (J2|, we have: 

1 



a — K K ( K m K d) / K d 

■ {K m - K d ) /K d + 1 J m m ($ ■ {K m - K d ) /K d + 1 

Eq. Q correctly predicts a small thermal conductivity of n m+ d ~ at q = 1, /9 « 0, and 
^ = (i.e., the cross-section is completely fractured), and a large thermal conductivity of 
K m+ d ~ K m at a ~ 0, /3 ~ and K d 7^ 0. 

Eq. Q can be simplified for common scenario where defect distribution is uniform and 
defect concentration is low. In such cases, K d 7^ (note that the K d = assumption used 
above implies a complete fracture of a cross-section normal to the heat flux direction) . When 
large d x , d y and d z values are used to contain the defects, we can always approach the limits 
K d ~ K m and (3 ■ (« m — Kd) / K d ~ 0. Using the relation 1/ (1 + x) « 1 — x for small x, Eq. 
Q can be approximated as 

K m+d & K m - (K m - K d ) ■ a ■ (3 = K m [l - (Km - K d ) / K m ■ p] = K m (1 - 7] ■ £) (5) 

where rj = f~ x ■ (K m — K d ) / K m is a constant for a given type of defect with a given size 
and given pattern (e.g., array). Eq. ^ indicates that thermal conductivity is a linear 
function of defect volume scattering density p. This verifies that for independent uniform 
defect distributions, distinction between areal and axial defect densities is not necessary. 
Furthermore, Eq. (J5| also indicates that thermal conductivity is a linear function of defect 
volume fraction £, for a given defect size, by way of the volume ratio /, Eq. ([I]). Eq. ^ 
can describe the scaling for the defect distribution shown in Fig. [IJa) if defect density is 
changed by the defect spacing 5. It can also be used to describe the scaling for the defect 
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distributions shown in Figs. [T^b) and Fig. [IJc) if the defect arrays remain unchanged (i.e., 
l x , L y , l z are kept constant) and the defect density is only changed by the spacing 5 between 
arrays. In the latter case, we are studying the scaling of arrays rather than the scaling of 
individual defect. 

Eq. ([5]) can also be derived based upon Fig. [6|c) as is described in Appendix [B| 

IV. MOLECULAR DYNAMICS VERIFICATION 

In this section, we verify the model by performing large scale MD simulations (400 proces- 
sors or more) using the bulk configuration (periodic boundary condition in the y— direction). 
Two series of MD simulations are conducted to explore effects of the areal and axial defect 
densities respectively. For areal effects, we use the configuration shown in Fig. |T|(b) where 
defect spacings in the x— and z— directions are small but are kept fixed whereas spacing 
between the x — z defect arrays in the y— direction is large and varied to change defect 
areal density. For the axial effects, we use the configuration shown in Fig. [I|c) where defect 
spacings in the y— and z— directions are small but are kept fixed whereas spacing between 
the y — z defect arrays in the x— direction is large and varied to change defect axial density. 
Because only the spacing in the direction with a sparse defect distribution is varied, this 
study satisfies our model assumption and can be used to test the model. The relatively small 
defect spacings in the other two directions enable small systems to be used to significantly 
reduce the computational expenses. In particular, all of our simulations used a small fixed 
dimension of n.3 = 6 (W ~ 19A) in the z— direction. 

A. Effects of Areal Defect Density 

To explore the areal defect density effect, our first series of simulations employ a fixed 
x— dimension (aligned with the flow of heat) of n\ = 136 (2L ~ 707A), and various y— 
dimensions between n 2 = 30 — 100 (t ~ 166 — 553A). As a reference of the dimensions, the 
smallest and the largest systems contain respectively 195,840 and 652,800 atoms. A heat 
flux of 0.00035 eV/A 2 ■ ps is used to introduce the temperature gradient. The choice of 
heat flux used in the work is to ensure that temperature difference between heat source and 
heat sink does not exceed 10 K but is significant enough to enable converged calculations. 
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Four equal-spaced voids are introduced between heat source and sink (the system contains 
a pair of heat sources and sinks and therefore eight voids). Notice that the void spacing 
in the x— direction does not affect the analysis of defect density in the y— direction, nor 
does the particular choice of four voids affect the generality of the results. As has been well 
established in the past [331 03], the use of different y— dimensions above a certain threshold 
does not affect the thermal conductivity of a defect-free matrix with periodic boundary 
conditions in the lateral directions. It does affect the distance between periodic defects in 
the y— direction (i.e., the spacing between the x — z defect arrays) and hence the areal defect 
density. Clearly, the defect spacings 6 X , 5 y , and 5 Z relate to system dimensions as 5 X = L/4, 
5 y = t, and 5 Z = W. Since L and W are not varied (i.e., we are not exploring the scaling 
in the x— and z— directions), we can simply set the dimensions of the (four) scattering 
volumes in these directions to their maxima d x = L/4 and d z = W. A large d y value is 
desired to contain the defect scattering region in the y— direction, but the geometry requires 
that d y < S y . We can maximize the constant d y by setting d y = 5 y , min = t min , where t min 
is the minimum thickness of all the computational systems used in the series. This ensures 
that the geometry constraint is satisfied for all the samples. Under these conditions, the 
defect densities a = t m i n /t, (3 — 1, and p = t min /t. Note that eight voids correspond to the 
removal of 96 atoms. For the smallest system, p — 1, and the defect volume fraction (site 
fraction) f = 96/195840 « 4.9 x 10" 4 = 4.9 x 1(T 2 %. This means that the defect volume 
fraction conversion factor for the series of samples is / = £/p = 4.9 x 10~ 2 %. 

According to Eq. ^ or (B5), thermal conductivity is a linear function of p. This 
relationship can be verified using data from direct simulations which are shown in Fig. [7] 
as a function of defect scattering density p or defect volume fraction £ using unfilled circles. 
It should be noted that the shaded area corresponds to the low defect density regime that 
cannot be easily calculated via MD using current computers. The thermal conductivity at 
p = 0, however, can be estimated using a smaller, defect-free system with periodic boundary 
conditions. The value thus obtained is shown in Fig. [7] with an unfilled star to distinguish 
it from other data points. The solid line in Fig. [7] is a linear function fitted to the simulated 
data. Note that the dashed line and unfilled diamonds show the effects of the axial defect 
density obtained in another series of MD simulations, to be discussed in the following section. 
At this point, the most important result in Fig. [7] is that the solid trend line fitted to the 
defect data predicts closely the defect-free value. This is a strong verification of Eq. ^ or 
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FIG. 7: Thermal conductivity as a function of the defect density p or £. Unfilled circles and solid 
line correspond to planar void arrays spaced 88 x 19 A 2 in the x — z planes, and unfilled diamonds 
and dashed line correspond to planar void arrays spaced 110 x 19 A 2 in the y — z planes. Spacing 
between arrays is varied to change defect densities. Note that the unfilled diamonds are not the 



direct result of MD simulations (see Section IV B) and hence are not associated with error bars. 



(B5). Also, since data at £ = and large £ can both be directly calculated using MD, Eq. 



(§ or (B5J can be used in an interpolation to reliably predict thermal conductivities in the 
low defect density regime that is not directly accessible by MD simulations. This is unlike 
the extrapolation typically used to infer the bulk limit in the presence of size-effects |33j . 

Numerically, our calculations indicated that the defect-free thermal conductivity at a 
sample length of n x = 136 (2L « 707A) is K m ^ 54 = 38.56 W/K ■ m (see Fig. Q. Note that 
thermal conductivity depends on both defects and sample length [331 1121 SSI and 
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as a result, we use K m ,354 instead of K m to indicate that the matrix thermal conductivity 
determined here pertains specifically to a sample length of L w 354 A. The void configura- 
tions discussed here are essentially planar void arrays spaced 88 x 19 A 2 in the x — z planes 
with spacing between arrays being varied between 166 and 553 A. For such particular defect 
configurations, we obtained a (non-dimensional) scaling coefficient of rj = 244.9 for the effect 
of defect volume fraction on thermal conductivity. 



B. Effects of Axial Defect Density 

To explore the axial defect density effect, our second simulation series employs a fixed 
y— dimension (perpendicular to the flow of heat) of = 20 (t llOA), and various x— 
dimension between n\ = 152 and 500 (2L ~ 790 — 2600A). As a reference of the dimensions, 
the smallest and the largest systems contain respectively 145,920 and 480,000 atoms. A heat 
flux of 0.0002eVyA 2 • ps is used to introduce the temperature gradient. One void is created 
in the middle between the heat source and sink (two voids in the system). Unlike the areal 
case where the change of thermal conductivity due to the change of system thickness t in the 
y— dimension comes only from the change of the areal defect density a (scales with 1/t), the 
change of thermal conductivity due to the change of system length L in the x— dimension 
comes from both: (a) the change of the axial defect density (5 (scales with 1/L) and, (b) 
the size-dependence of the interfacial scattering [33], E2J H3j I54T - I56] . Our scaling model can 
be used to derive an analytical expression of thermal conductivity as a function of both 



defect density and sample length, as described by Eq. (CI) in Appendix C Based upon the 



relation p oc 1/L, Eq. (CI) indicates that thermal resistivity of defective material with a 
finite length L is a linear function of 1/L. The validity of the analytical model on the axial 
defect density can therefore be verified by checking this linear relationship. Furthermore, 



since the length scaling coefficient p (refer to Eq. CI) can be determined from a series of 
defect-free samples with different lengths [33], the effect of defects can be isolated. 

Here, the simulated thermal resistivity results are shown in Fig. [8] as a function of 1/L 
using the unfilled diamonds. For comparison, similar data for defect-free samples obtained 
previously [33] are included as the filled circles. Note that the lines are calculated from a 



function fitted to Eq. (CI) as will be described in the following. It can be seen that the 



linear relations are very well satisfied. Most importantly, the thermal conductivity data for 
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FIG. 8: Thermal resistivity as a function of inverse of sample length 1/L with and without the 
defect in the middle of sample. 

the defective samples is seen to be lower than that of the defect-free samples, and both 
defective and defect-free samples approach the same K mj0 o limit at L — > oo (also p — > 0). 
As described in Appendix [AJ a direct test of how good the linear relationship extends to 
the infinite sample length is not possible using merely defect-free samples because it is 
increasingly difficult to obtain accurate thermal conductivities from MD simulations when 
sample length is increased. However, the convergence of defective and defect-free samples 
to the same thermal conductivity at the infinite sample length limit is a strong verification 



of the linear scaling law, Eq. (CI). Eq. (CI) therefore can accurately predict thermal 



conductivities in the defect density and sample dimension space that cannot be directly 
assessable by MD. 

Now, we compare the results of this section to those of the previous section which ex- 
amined the influence of areal density. Here, we set d y = 5 y = t, d z = S z = W. In the 
previous section, a defect volume scattering density p — 1 corresponds to eight voids in a 
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136 x 30 x 6 cells 2 system. This density is equivalent to two voids in a 51 x 20 x 6 cells 2 
system. To match the previous defect volume scattering density definition, we hence choose 
d x = 51/2 cells = 132.6 A (< 5 Xjmin = L mini where L min corresponds to the minimum x— 
dimension per defect used in the simulation series). At these values, the defect densities 



a — 1, /3 — 132. 6/L, and p = 132. Q/L. Using these definitions, we fit Eq. (CI) to our data. 
This reproduced the bulk limit of the thermal conductivity value of K mj0 o = 184.97 W/K -m 
obtained previously [33], and resulted in the determination of a length scaling coefficient of 
p = 7.4248 x 10 -10 K ■ m 2 /W, and a defect volume scattering density scaling coefficient of 



q = 0.0109 K ■ m/W, see Eq. (CI). The lines shown in Fig. [8] are calculated using these 
parameters. Note that in this series of simulations, the defects are essentially planar arrays 
of voids spaced 110 x 19 A 2 in the y — z planes with array spacing being varied between 395 
and 1300 A. 



Eq. (CI) also allows us to cast the thermal conductivity obtained at one sample length 



Li to another L 2 at a constant defect density: 

1 1 ( 1 1 \ 

r^ = rm+p* rr _ T ( 6 ) 

K m+d{P)\L 2 K m+d \P) \L 1 ^2 J 

Using Eq. ([6]), thermal conductivity vs. defect density data at a fixed sample length of 
Hi = 136 (2L ~ 707 A) is obtained. Fitting Eq. (pi) or (B5) to such data resulted in 



^m,345 = 38.56 W/K ■ m and a scattering coefficient of rj = 918.4 for the effect of defect 
volume fraction on thermal conductivity. Pertaining to the same sample length as the 
unfilled circles in Fig. [7j the converted data and the corresponding fitted function are shown 
as unfilled diamonds and dash line in Fig. [7j It can be seen that like the areal defect density 
effect, reducing the axial defect density also causes the thermal conductivity to approach 
the value of the defect-free sample. With consideration of Appendix [A] which discusses 
the difficulties in directly testing the thermal conductivity at very large sample length, the 
convergence to the same defect-free sample point from two defect relations shown in Fig. [7] 
is again a strong verification of the analytical model. 

For non-zero defect densities, thermal conductivities obtained from the areal and axial 
series are different, with different coefficients rj. This means that when defects are not sparse, 
their scattering regions can overlap, resulting in different effects on thermal conductivities 
depending on how defects are distributed. For example, the configurations shown in Figs. 
[TJb) and [T^c) have different scaling coefficients. 
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FIG. 9: Thermal conductivity of a thin film as a function of void-surface distance predicted by MD 
simulations and an analytical model. 

V. EFFECT OF DEFECT SPATIAL DISTRIBUTION 

A series of MD simulations are performed to study the effect of defect location with 
respect to a surface using a film configuration with the free boundary condition in the y— 
direction. The system has a fixed dimension of n\ = 136 (2L rs 707 A), ri2 = 20 (t ~ llOA), 
and n 3 = 6 (W ~ 19 A) for a total of 130,560 atoms. Four voids are created between the 
heat source and sink, resulting in a void volume fraction of £ = 7.35 x 10~ 2 %. These defects 
are uniformly distributed in the x— direction (aligned with the heat flow), but the position 
of the row of defects is varied along the y— direction so that it has different distance s from 
the free surface, Fig. [2] The calculated thermal conductivities as a function of s are shown 
with unfilled diamonds in Fig. [9] Fig. [9] clearly indicates that thermal conductivity depends 
on the defect-surface distance. In particular, thermal conductivity increases as defects move 
towards the surface. 
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FIG. 10: Illustration of defect population model. 



Analysis can be used to understand the results shown in Fig. [9j Eq. ^ or (B5) indicates 
that if the total scattering strength of one type of defects is T]\ (aggregating the dependence 
on the defect volume fraction £), then the thermal conductivity K\ is reduced from the matrix 
value kq through k± = Kq • (1 — r]\). Similarly, if there is a second type of defects with the 
scattering strength 772, then the thermal conductivity K2 of the material can be viewed as 
reduced from the new matrix value k% through K2 = «i • (1 — ^2) = ^0 ' (1 — Vi) ' (1 — ^2)- In 
general, we can assume a multiplication rule of k n = Ko-H^ (1 — rji) to account for iV types 
of defects. To relate to the MD results, a two dimensional illustration of the system is shown 



in Fig. 10 , where voids are assumed to lie on the center line of the shaded region a distance 



s below the surface, and the total sample thickness is t. Because the interaction between 
surface and defects needs to be addressed, the thermal conductivity must be considered 
locally as a function of position y. 

It is obvious from Fig. [3] that the local scattering strength r]^ s (d) of a surface is a 
decreasing function of distance d from the surface. We postulate that such a behavior can 
be well captured by the function: 



%s (d) = S K ,s ■ exp (-// ■ d) 



(7) 



where S KtS and [i are constants. The local thermal conductivity Ki >s (y) due to the presence 
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of the surfaces (but no defects) is therefore: 

(v) = «m • [1 - *7j,« (*/ 2 - 2/)] ■ [1 - Vi,s it/2 + y)\ (8) 

where t/2 ± y are the distances of a local site y from the two surfaces. Similarly, we can 
assume that the local scattering strength rji^ (d) of defects is a decreasing function of distance 
d from the defect plane based upon Fig. [3j Again, we postulate that the behavior can be 
well captured by the function: 

rji,d (d) = 6 Kid ■ exp {—v ■ d) (9) 

where and v are constants. The local thermal conductivity (y) due to the presence 
of the defects (but no surface) is then 

Kj, d (y) = Km • [1 - 77j lrf (|t/2 - s - y|)] (10) 
where \t/2 — s — y\ measures the distance of a local site y from the defect plane. Note that 



Eq. (10) prescribes the "apparent" thermal conductivity in a thin slice parallel to the defect 



plane at a location of y. For such a slice, Eq. (10) correctly specifies a drop of the thermal 
conductivity by a fraction of 5 K ^ at the defect plane y = t/2 — s, and the recover of the 
bulk value K m when y is far away from the defect plane with the parameter v essentially 
capturing the scattering distance. 

In regions where both surface and defects have significant scattering, the total local 
thermal conductivity K\ (y) can be written as: 

ki (y) = K m ■ [1 - m>a (t/2 - y)] ■ [1 - m , s (t/2 + y)] ■ [1 - m>d (\t/2 - s - y\)\ (11) 



The apparent thermal conductivity of the system can then be found by averaging Eq. (11) 

as 

i r y=t * 

k = - / Kl {y) dy (12) 



t 



y= 



2 



To apply Eq. (12), parameters K m , fi, u, S K>S , and are needed. For the simulated 
sample length, the matrix conductivity is n m = /t m ,354 = 38.56 W/K ■ m. The other param- 
eters can be estimated from independent simulations. First, we express the average thermal 
conductivity of a system containing the surfaces but not the defects as 

2 ry =t / 2 

K s = - I K i,s (y) dy (13) 

1 Jy=0 
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Here we only integrate half of the sample thickness due to the symmetry of the problem. 
A series of MD simulations are carried out to obtain n s vs. t data for defect-free thin film 
samples. In particular, we use a fixed x— dimension of n\ = 136 (2L « 707A) and various 
y— dimensions of 712 = 30, 50, 70, 96 with free y— surfaces (t ~ 166 — 530A). In addition, we 
perform an additional simulation with the periodic condition in the y— direction at n 2 = 10 



(t — > 00). By fitting Eq. (13) to the k s vs. t results obtained from MD simulations, we find 
S KjS = 0.15135 and n = 0.00705 A -1 . 

The average thermal conductivity of a periodic (i.e., no surfaces) system containing de- 
fects on the center line (i.e., s = t/2) is expressed as 



2 rv =t / 2 

K d = - (y) dy (14) 

ly=0 



t 



Note that the defect scattering densities used in Fig. [7] (unfilled circles) satisfy (3 = 1 and 
P — tmin/t, where t m i n = 166 A (n 2 = 30) is the minimum system thickness used in the 
series. Hence, the thermal conductivity vs. defect density (p) data shown in Fig. [7] can be 
converted to the thermal conductivity vs. \jt data. By fitting Eq. ( 14 ) to these conductivity 
vs. l/t data, we obtain 5 M = 0.29017 and v = 0.02759 A -1 . 



Based upon the parameters thus obtained, Eq. ( 12 ) is used to calculate thermal conduc- 
tivity as a function of s, and the results are included as the dots in Fig. |9j A good agreement 
between the analytical prediction and the MD data is clearly shown, thereby verifying the 



analytical model. The mechanism of the defect spatial effects is now clear. Eq. (11) indi- 
cates that when s — > 0, the defect affected region merges with one of the surface affected 
regions. Consequently, the total scattering affected region is reduced from the case where 
defect is in the middle of the sample (s — > t/2). It is this consolidation of the defect and 
surface scattering that causes an increase in thermal conductivity when defects approach 
the surface. 



VI. DISCUSSION 

Thermal conductivity as a function of void density, size, and distribution have not been 
experimentally measured. However, analytical expressions for thermal conductivity as a 
function of porosity (i.e., volume fraction £) have been experimentally derived for a variety 
of porous materials [20TI24"] . While different forms of analytical expressions are used to 
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enable good fit of experimental data up to a high porosity of 0.8 [23], all expressions can 
be reduced accurately to a linear function n m+ d = K m - (1 — r\ ■ £) when porosity £ is small 
enough (below 0.05). This is in good agreement with Eq. ^ or (B5), which is also valid for 
void volume fraction far less than 0.05. However, there is a significant difference between 
the nanoscale void effects and the macro-scale porosity effects as the scaling parameter r] for 
different porous materials falls between 1.0 and 4.6 [20TI24"] whereas 77 is 244.9 and 918.4 for 
the two void configurations explored in Fig. [7} This means that the thermal conductivity 
at a given void density cannot be interpolated linearly between a defect-free sample and 
porous materials, i.e. as void size changes. 

The discussion of the spatial effects of defects presented in the preceding section indicates 
that thermal conductivity increases when defects move closer because their scattering regions 
overlap resulting in a reduction of total scattering volume. This can account for the difference 
between voids and pores, which can be more clearly illustrated using Eq. ([I]). When a large 
number of voids are closely packed to form a pore, the pore size £ x , £ y , and £ z become 
very large (> 2/im [20]) but the scattering dimension £ s remains small (assumed to be 
comparable to that in the void case). This means that the defect volume-to-scattering 
density conversion factor / approaches 1 for large pores. If the thermal conductivity i^d of 
a pore approaches 0, rj approaches 1 by definition. This explains why the parameter 7] for 
a variety of porous materials falls in the order of 1. On the other hand, if a large pore is 
split into a large number of voids of sizes around 5 A distributed uniformly in the material, 
a large number of independent scattering volumes will be created, resulting in a significant 
reduction of the conversion factor /, a significant increase in the parameter 7], and hence 
a significant reduction of thermal conductivity. This accounts well with the previous MD 
results that vacancies cause a more significant reduction of thermal conductivity than voids 
given the same defect site fraction [38]. Interestingly, experiments also indicate that at a 
given porosity, reducing pore sizes causes an increase in rj [20], which is consistent with an 
increase in total scattering volume. Our analytical model and MD data, therefore, are well 
corroborated by the experimental data on porous materials. 

At the lower limit of defect size, point defects have been studied extensively and the 
low temperature effect of Rayleigh scattering is well-known. However, there have been 
relatively few experimental studies of actual vacancies - most experimental data is for isotopic 
substitutions. Che et al. [57] calculated the effects of point vacancy concentrations in the 
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range 0.01-0.16% for carbon. Unfortunately, the power-law-like empirical model they fitted 
to the data does not have finite derivative at a zero defect concentration. Nevertheless, 
it is clear that this limit of phonon scattering from voids results in a scaling coefficient 
significantly larger than the macroscopic porosity case. 

Molecular dynamics simulations have been recently used to study thermal transport of 
nanoporous crystalline [37] and amorphous [36] silicon. The results clearly indicate that 
thermal conductivity depends on pore size and pore fraction. In particular, thermal conduc- 
tivity was found to reduce with an increasing interfacial area concentration, which relates 
well to the effects of the scattering volume discussed above. Our studies, therefore, verify 
the previous results. In addition to the pore size and pore fraction effects, we further show 
that thermal conductivity is sensitive to defect population (for example, areal and axial 
populations have different effects). 

VII. CONCLUSIONS 

An integrated approach combining a physically-motivated analytical model, large scale 
MD simulations, and extensive experimental thermal conductivity data of porous materials 
is used to study defect effects on thermal conductivity. Corroborated results lead to an 
explicit functional expression of thermal conductivity on defect density, size, and spatial 
population. The following conclusions are of particular interests to both theoretical thermal 
transport studies and phonon engineering of materials: 

(a) Thermal conductivity depends strongly on total scattering volume of defects. This 
scattering volume differs from the physical volume of defects. It, however, can be 
linearly correlated with physical volume of defects when defect type and size are fixed. 
It is the defect configuration that minimizes this scattering volume but not the physical 
volume that will increase thermal conductivity. 

(b) When defects are close, their scattering regions overlap, resulting in reduced total scat- 
tering volume. As a result, thermal conductivity increases when voids move towards 
surfaces, or small voids collapse to form large pores. 

(c) For uniform, sparse defect distribution with a given defect size, thermal conductivity is 
a linear function of defect volume fraction. However, thermal conductivity at a given 
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void density cannot be interpolated between defect-free samples and macroscopically 
porous materials due to the difference in defect sizes. In general, the dependence of 
thermal conductivity on defects is not simply through volume fraction, but strongly 
depends on the size distribution of the voids. 

(d) The analytical model enables thermal conductivities obtained from molecular dy- 
namics simulations to be extrapolated/interpolated reliably to realistic defect density 
ranges, as well as the defect-free limit. 

Finally, we point out that the success of our approach can be related to the explicit 
incorporation of scattering of short wavelength and short mean free path phonons in direct 
molecular simulations. These modes are scattered the most given the general behavior of 
Rayleigh scattering and therefore their contribution to the overall conductivity is the most 
sensitive to changes in defect density. On the other hand, the longer wavelength modes are 
less affected by point-like defects and their contribution to the thermal conductivity is given 
by the scaling analysis that estimates the long sample length limit. 
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Appendix A: Examination of the Length Scaling Law 

Michalski [54J first discovered that when a heat flux passes through a material with a 
length L, the inverse of thermal conductivity of the material, 1/k, linearly increases with 
1/L (i.e., 1/k — l//«oo + p/L where is the thermal conductivity at L — > 00 limit and 
p is a length independent slope parameter). This scaling law was later also reported by 
other researchers [SSI EH]- Since then, such a law has been used to calculate by linearly 
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extrapolating the thermal conductivity data obtained in "direct method" MD simulations at 
a series of short sample lengths [32, 33, 52J . However, controversy has arisen in recent studies 
regarding the validity of the linear extrapolation method. In fact, our MD simulations 
on GaN thermal transport using periodic boundary condition in the heat flux direction 
[33] indicated that while the linear relation holds well at low temperatures of 300 K and 
500 K, it deviates at 800 K. In particular at 800 K , the slope increases discontinuously 
when the sample length is increased above a threshold value, indicating two regimes and 
resulting in a larger extrapolated thermal conductivity when the sample length range is 
above this threshold. Similar phenomenon was later observed by other researchers [58- 
[61] . Based on the studies of Si using non-periodic boundary conditions in the thermal 
transport direction and motivated by the concept that diffusive thermal transport does not 
hold at length-scales below the phonon mean free path, Sellen et al [58] suggested that the 
linear extrapolation procedure might not be valid unless the sample lengths used in MD 
simulations are significantly longer than those commonly used in literature. Recently, we 
derived an analytical expression relating thermal conductivity to dimensions in all three 
coordinate directions, and verified using large scale molecular dynamics simulations that 
this expression is satisfied when the sample dimensions are above certain threshold values 

32JH3]. 

While the literature studies cited above involve some of the largest MD simulations per- 
formed so far in the field, they have not given a fully consistent story. For example, Sellen 
et al's work [58] suggested that the extrapolated bulk thermal conductivity n^, is underes- 
timated if the MD systems are too short, and that realistic /too can eventually be achieved 
if the MD samples are sufficiently long. However, our previous work [33] indicated that at 
least for GaN at 800 K, unrealistically large be obtained using MD sample lengths 

in a range shorter than that proposed by Sellen et al's work [58] and reasonable 
obtained only when sample length is further reduced. While the underlying mechanism for 
this abnormal phenomenon was not clear, we emphasize that our results were obtained from 
very long averaging times (at least 40 ns and some are significantly longer) and hence the 
unrealistically large extrapolated thermal conductivity was unlikely to be caused by statis- 
tical errors. On the other hand, Yang et al [59] indicated that for Si with sample lengths 
up to ~ 200 nm, the linear relation holds well under the periodic boundary condition in 
the thermal transport direction, but changes discontinuously under a non-periodic condi- 
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tion. Clearly, the sample length range within which the linear relation holds depends on the 
boundary conditions. To help elucidate the literature observations on the linear scaling law, 
here we re-examine effects of both the boundary condition and sample length by perform- 
ing direct method MD simulations using averaging times and a sample length significantly 
beyond those used previously The same approach used previously [33] is followed to study 
the perfect GaN system at a 300 K temperature. 

First, simulations are performed to calculate thermal conductivities of GaN crystals at 
different sample lengths. The crystal configurations are the same as those used previously 
[55] except that free boundary conditions are used in the thermal transport direction here 
whereas the previous work employed a periodic boundary condition. Two scenarios are 
considered, one with the two free-ends fixed, the other with the two free-ends free. To 
ensure highly-converged results, we use a very long averaging time of 110 ns. The results 
of thermal resistivity (1/k) vs. inverse of thermal transport distance (1/L) obtained for 



the fixed and free ends are shown respectively as blue circles and red squares in Fig. 11 



For comparison, the data obtained previously for the periodic boundary conditions [33J are 



included in Fig. 11 using unfilled diamonds. Fig. 11 clearly shows that the boundary 
condition has a significant effect on the linearity given the same sample length range. For 
the GaN system explored here, the free-end non-periodic condition has the most significant 
non-linearity, whereas the periodic condition produces the most linear relation. 

One obvious question is how the phonon mean free paths of the collection of phonons 
affect the linear scaling law. The phonon mean free paths (and their associated relaxation 
times) are very difficult to calculate from MD simulations [531 162] , but the mean free path 
for GaN at 300 K has been experimentally measured to be around 0.1 \im [16]. Hence, a 
GaN crystal with a sample length of 2L = 1.3 fim is used to perform an additional MD 
simulation under the periodic boundary condition. While the thermal transport distance 
L equals half of the sample length under the periodic boundary condition and the phonon 
mean free path is uncertain under the simulated condition, it should be noted that our L 
value is comfortably above the experimental phonon mean free path. 

An increased sample length causes a significant increase in both the equilibration time for 
establishing temperature gradient and the averaging time for highly-converged results (at 
least scales with L 2 ). We hence perform a very long pre-MD run (64 ns) to equilibrate the 
temperature gradient, and use a very long (32 ns) simulation over which average tempera- 
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FIG. 11: Thermal resistivity of GaN crystals as a function of inverse of heat source-sink distance 
(1/L) at 300 K temperature. 

ture profiles are computed and used to evaluate the thermal conductivity. The result thus 



obtained is included in Fig. 11 using the filled circle. It can be seen that the filled circle is 
very close to the solid line deduced from a linear expression of the MD data simulated under 
the periodic boundary condition in a small sample length range (i.e., the unfilled diamonds). 
Fig. [TTJ therefore, does not indicate that the length scaling law is violated for GaN crystals 
under periodic boundary conditions. This is a result that corroborates the findings in Yang 
et al [32]. However, both our previous GaN work at 800 K |33j and the calculations here 
strongly indicate that converged thermal conductivity at a long sample length is extremely 
difficult to obtain. Note that the simulation time we use here is limited by the computing 
resources currently available to us and we only calculate one point at the very long sample 
length. While our MD time is significantly longer than that of most work reported in liter- 
ature, we propose to significantly increase the averaging time in the future to re-calculate 
many points in a long sample length range when more powerful computing resources are 
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available. 



Appendix B: Second Derivation of Thermal Conductivity as a Function of Defect 
Density, Size, and Distribution 



First, the material volume 5 X ■ 5 y ■ 5 Z is divided along the x direction into three slabs, the 
central slab has a thickness d x that contains the defect, and the remaining two end slabs 
have a total thickness of S x — d x , as shown in the left frame of Fig. |6](c). The apparent 
thermal conductivity of the central slab is assumed to be K dx , and the thermal conductivity 
of the remaining two end slabs remains to be the matrix value K m . Because the heat 
conducts through the three slabs in serial, the overall thermal resistivity of the material can 
be represented as a length weighted average: 



$x K d x $x 



K, 



K d x K T 



(Bl) 



In Eq. (Bl ), re^ can be further expressed. This is illustrated in the right frame of Fig. |6[c), 
where the central slab d x ■ 5 y ■ 5 Z is decomposed into a central volume d x ■ d y ■ d z and the 
remaining material. Because heat conducts through the central and remaining materials in 
parallel, the overall thermal conductivity of the central slab is an area weighted average: 



Kdx 



dy ' d z dy ' 6 Z dy • d z 

dy ■ d z 



5 V ■ 5 Z 



a ■ Kd + (1 — a) ■ K r 



Substituting Eq. (|B2|) into Eq. (|B1|), we have: 
1 1 ] 



^m+d 



K, 



1 



1 



K, 



1 



( K d K m) / K r 



a - (K d - K m ) /n m + 1 



(B2) 



■P (B3) 



a ■ (K d - K m ) /K m + 1 

Like Eq. Q, Eq. (B3) also correctly predicts a small thermal conductivity of K m+( i « at 
a = 1, (3 ~ 0, and Kd = 0, and a large thermal conductivity of K m+ d ~ K m at a ~ 0, (3 « 0, 
and Kd 7^ 0. Again for a low density of uniform defects, we can choose large d x , d y and d z 
values so that K d ~ K m and a ■ (Kd — K m ) / K m « 0. Eq. (B3) can then be simplified as 

1 



^m+d 



Eq. (B4) is further written as 

t 

^rn+d 



1 P ' (^d K m) / 

Eq. (B5) is the same as Eq. 



(K m - K d ) ■ p = K m (1 - rj ■ £) 



(B4) 



(B5) 
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FIG. 12: Effects of defect and end surfaces on thermal transport. 

Appendix C: Combined Effects of Defects and Finite Length between Heat Source 
and Sink-surface 



In the axial defect density series of MD simulations, samples of different lengths L were 
used to change the defect density. However, the resulting thermal conductivity is affected 
by both sample length [33J H21 ffiH and defect density. Here we derive a thermal 

conductivity expression to separate the length and defect effects. 



As shown in Fig. [T2j we assume that heat flows through a medium with a finite length 
L that contains a defect in the middle. If we use a large d x to fully contain the defect effect 
and a large u to contain the end surface effect, then the overall thermal conductivity of the 
defective section and the end surface regions can be taken as constants ft^ and ft e respectively 
and the conductivity of the remaining material is that of the infinite bulk ft mj00 . Note that 
because both defect and end surface effects are considered here, we use ft m)00 instead of ft m 
to emphasize that the thermal conductivity of the matrix is free of both defects and end 
surfaces (i.e., the length L is infinite). Because heat flows through different regions in serial, 
the overall thermal resistivity is a length weighted average: 

1 2w 1 (LI L-2u-d n . 1 



l 4 

L k p L 



1 



2ljJ (ft mj0 o K e) d x (ft 
L ' ftm. rm ' ftp 



1 p 

+ y + q-p 



L ■ ft r 
1 



+ 



m,oo 
1 



p 



(Cl) 



where p 



2(jJ ■ (ft m ,oo K e 



ft, 



• fte, q 



ft 



in. to ^d) ' K m ,oo ' K d 



1 are constants for 



given types of surfaces (or interfaces) and defects, and p is taken to be p = d x /L assuming 
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a = 1 (i.e., we can use the entire cross-section area to contain the areal scattering effects of 
defects). Eq. (CI) can be used to separate defect and system length effects. 
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